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Abstract 

We study the ergodic properties of a two-dimensional self-gravitating system using molecular 
dynamics simulations. We apply three different tests for ergodicity; a direct method comparing 
the time average of a particle momentum and position to the respective ensemble average, sojourn 
times statistics and the dynamical functional method. For comparison purposes they are also 
applied to a short-range interacting system and to the Hamiltonian mean-field model. Our results 
show that a two-dimensional self-gravitating system takes a very long time to establish ergodicity. 
If a Kac factor is used in the potential energy, such that the total energy is extensive, then this 
time is independent of particle number, and diverges with '/N without a Kac factor. 
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I. INTRODUCTION 


A^-body systems with long-range interactions have some peculiarities with respect to 
systems with short-range interactions, having drawn much attention at least along the last 


two decades HHsi. and for a longer time if one considers self-gravitating systems and charged 


plasmas. Starting from a non-equilibrium conhguration, a short-range interacting system 
evolves to thermodynamic equilibrium in a relatively small relaxation time, while a long- 
range interacting system evolves over different stages characterized by time-scales differing 
by orders of magnitude, and taking a very long time to reach thermodynamic equilibrium 
for a hnite number of particles. The initial stage of evolution is a violent relaxation into a 
Quasi-Stationary State (QSS) in a time roughly independent on the number of particles [6]. 
The relaxation time to thermodynamic equilibrium diverges in the N ^ oo limit (Vlasov 
limit) and gets trapped in a non-equilibrium non-Gaussian stationary state. 

A pair interaction potential is long-ranged if it decays at long distances as r~'^ with 7 < d, 
with d the spatial dimension and r the inter-particle distance. This implies that all particles, 
no matter how far, contribute to the total energy. Consequently the system is non-additive, 
violating the fundamental axiom of thermodynamics: non-additivity of entropy. It is worth 
remembering that this is not in contradiction with the second law of Thermodynamics, which 
is always valid . Another common consequence is non-ergodicity, most extensively studied 
for the Hamiltonian mean held model [HI nHn], but also for one and three-dimensional 
self gravitating systems [TUl - 125] . Quite curiously, and up to the authors knowledge, no 
previous detailed study was devoted to the ergodicity of the two-dimensional case, at least 
not in the sense discussed here. In the Vlasov limit, the dynamics of the system is exactly 
described by the Vlasov equation, which is essentially the Liouville equation for the one- 
particle distribution function evolving in the mean-held due to all other particles, and given 


by [21I2H]: 



( 1 ) 


where /(p, r; t) is the one-particle distribution function, p and r the momentum and position 
vectors in d spatial dimensions, respectively, and 



( 2 ) 


is the mean held force at position r and time t with V(r — r') the inter-particle potential. For 


2 


finite N the right-hand side of Eq. Q does no longer vanish and is given by the pertinent 
collisional contributions (granularity effects) |2S]- For a long-range interacting system the 
thermodynamic limit is ill dehned [2Z], and the Vlasov limit can be taken consistently by 
introducing a Kac factor 1/N in the interaction potential, such that the energy becomes 
extensive, but still non-additive, and that the iV —>■ cx) limit converges |28j . 

The study of ergodicity was pioneered by Boltzmann in his works on the foundations 
on Statistical Mechanics [29], and latter extended by Birkhoff [30] and Khinchin [31]. The 
ergodic hypothesis states that time average of a dynamical functions b{x) equals the ensemble 
average: 

lim ^ ^ dt b{xt) = b{x) = {b{x)) = /dpo b{x), (3) 

tf^oo tf — to Jto J 

where x denotes a point in the system state space (or phase space for a Hamiltonian classical 
system), b{x) is the time average of &(x), {b{x)) its ensemble average and dpo a statistical 
measure. Proving ergodicity rigorously is a difficult task and has been accomplished only 
for a few cases. Most studies rely on different methods such as determining the existence 
of gaps in phase-space mm, direct comparison of time and ensemble averages for the 
momentum variable [13], sojourn time statistics for cells in phase-space [TH] 1321 ES], testing 
for equipartition of energy [T^ 120] 131] and the dynamical functional approach [33H3Z]- The 
latter can also be used to determined weather the system is mixing, which is a stronger 
property than ergodicity. The dynamics of a system is mixing in its phase-space S if for an 
invariant measure p, T a map preserving fi and si,S 2 C S, we have that [39]: 

h(5i) /i(s2nT^Si) 

™ ^-• (4) 

H{S) k^oo /i(s2) 

Non-ergodicity can also be classihed as strong if some regions of phase-space are non- 
accessible to the system and weak if all regions are accessible but not equally visited. 

Ergodicity can be studied in two complementary ways. The hrst possibility is to consider 
the evolution in phase space of the V-particle system, with time averages computed over the 
trajectory of the point representing the state of the whole system. Ensemble averages are 
then obtained by considering all points in the energy hypersurface as equally probable, i. e. 
for the microcanonical equilibrium. This is equivalent to state that the total time spent in 

a phase-space region is asymptotically proportional to its measure. The other possibility is 

to consider single particle histories as realizations of a stochastic process. In this case time 
averages are taken over one particle history, while ensemble averages are taken over the set 
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of all N particles. This approach is particularly suited for a direct experimental verification 
even for short-range interacting systems [3H]- In this case ergodicity, i. e. the coincidence 
of the time and ensemble averages is not equivalent to thermodynamic equilibrium as will 
become evident below. We consider the latter case in the present work, with S in Eq. (|4]) the 
one-particle phase-space, for a system composed by a single particle evolving in a stationary 
potential, either in the thermodynamic equilibrium or in a (quasi-) stationary state. 

In this paper we hrst briefly study a gas of elastic hard-discs as an example of short-range 
interacting system, and then revisit the ergodic properties of the Hamiltonian Mean Field 
(HMF) model extending the results of Ref. [TB] to include non-homogeneous quasi-stationary 
states. Both system are used as comparison standards for our main interest here, a two- 
dimensional self-gravitating Wbody system. The HMF model is formed by N particles on 
a circle with Hamiltonian [8]: 

N 2 1 JV 

^ = + E (5) 

k=l ^ k,l=l 

where 6 k is the position angle on particle k and pk its conjugate momentum. It is a solvable 
model at equilibrium |H1 [9] at its dynamics can be simulated with a smaller computational 
cost if compared to other long-range interacting systems [10] and widely studied in the 
literature (see HHH and references therein). A two-dimension self-gravitating system with 
N identical particles, as described by the Hamiltonian: 

N 2 1 V 

= + ^ E >“(lr‘-n| + £). (6) 

k=l k,l=l 

where the logarithmic potential is the solution of the Poisson equation in two dimensions [12] , 
all masses are set to unity and e is a softening small parameter, commonly used in simulations 
of self-gravitating systems to avoid the divergence in the potential at zero distance IS]. Note 
that we use a Kac factor 1/A^ in the potential. 

The potential in Eq. ([^ is not upper bound and is therefore conhning, consequently 
avoiding the difficulty of particle evaporation in three-dimensional gravity. This model 
was used for instance in the study of anomalous diffusion in a collapsing phase [13], the 
determination of thermodynamic equilibrium properties m, collisional relaxation [25105], 
violent relaxation [B2l 06] and cooling in self-gravitating accretion discs 071 . 

The structure of the paper is as follows: In section [TT| we succinctly present the methods 
for testing ergodicity used in the present work, and apply then to a two-dimensional hard- 
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disc gas in Section and to the HMF model in Section W The ergodic properties of a 
two-dimensional self-gravitating system are discussed in Section |V| We close the paper with 


some concluding remarks in Section VI 


II. TESTING ERGODICITY 


Direct method 


We consider the time evolution of a single particle to compute time averages while ensem¬ 
ble averages are obtained by taking the average over the N particles in the whole system. 
For an ergodic system both quantities coincide. Let us then consider the momentum Pk{t) 
of particle k and its position Xk{t) at time t and write their time averages as: 


and 


Pkit) = -£pA:(jAt), 

^ i=0 


1 " 

^k{t) = - 

^ 3=0 


(7) 


( 8 ) 


respectively, where At is a constant time interval, that we take as being the numeric in¬ 
tegration time step and m is the total number of such intervals. The time averages are 
computed up to a given time t denoted by the argument of Pk(t) and Xk(t). On the other 
hand, ensemble averages at time t are given by 


N 


{pit)) = ^Y.Piit)^ 

i=0 


and 


1 


N 


{x{t)) = — 


(9) 


( 10 ) 


i=0 


For an ergodic system Pf, and Xk are the same for all particles, and therefore the standard 
deviations: 

CTp(t) = 


1 


N 




and 


(Jxit) = 


1 


N 




( 11 ) 


( 12 ) 


must vanish asymptotically with time. In Eqs. (0 and (0 we used {p{t)) = {1/N) EkPkit) 
and {x{t)) = il/N)J2kXkit)- 
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The direct method then consists in verifying numerically that dp and ax{t) do indeed 
vanish (or approach zero) for some sufficiently long time t [T3]. For more than one spatial 
dimension this must be true for each momentum and position components. 

It is interesting to note the connection of ergodicity breaking and anomalous diffusion in 
long-range interacting systems. In the limit At —)• 0 Eq. ([^ is rewritten as: 

1 /■* _ , , m r^kit) 


Pk = - dtpk{t) = — dxk = — [Xk{t) - a;fc(0)J. 
t Jo t JxM) t 


(13) 


Plugging Eq. (13) into Eq. 10 and assuming (p(t)) = 0 we obtain: 


CTp(t) = 


m 


N 


t\Nt, 


[Xkit) -Xfc(O)]' 


(14) 


Anomalous diffusion is ubiquitous in systems with long-range interactions, with the standard 
deviation of the distance traveled by a particle satisfying: 


1 ^ 

-Xk{o)f = cr, 

k=l 


(16) 


with C and p constants, and the case p = 1 corresponding to normal diffusion. From 
Eqs. (0, @ and ( [M] ) we have 


crp(t) = 


(16) 


As a consequence the system is always ergodic in the momentum variable whenever p <2. 
Nevertheless the observation time required to reach ergodicity can be very long for p smaller 
than but close to 2. This time also diverges if transport is ballistic or if C —)■ cx) for N —)■ oo. 
A similar approach is not possible for the position variable as in this case there is no relation 


analogous to Eq. (13). Anomalous diffusion in the HMF model was thoroughly studied 
in the literature 


], although some points still deserve more investigation. For 
self-gravitating systems see Ref. |M| and references therein. 


B. Sojourn time statistics 

This approach was introduced by Rebenshtok and Barkai for a Continuous Time Random 
Walk (CTRW) [32]. Consider a system of discrete states labeled by an integer index, each 
one visited intermittently for a given time. The time spent in state k during the j-th 
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(s) 

visitation is denoted by j and called a sojonrn time. The snm over all sojonrn times for a 
given state yields its residence time: 


tk = ^t 


(s) 

kj' 


j 


(17) 


The probability density fnnction for a weak non-ergodic system for the time average O of 
an observable O can be written as [32]: 


/(«)(0) 


— lim Im 

TT 


. eL,p?(o-Oi + «)”''^ J’ 


(18) 


where a is a constant in the interval (0,2], Ok is the value of the observable in state fc, 
is the probability for the system to be in the state k and L the number of discrete states. 
The time average of the observable O is thus: 


0 = ^j2tk0k, (19) 

hot k=i 


with hot = Yli=k ^k- If the statistics of the sojourn times have a power law tail, then the 


system is non-ergodic. Otherwise the system is ergodic and we have that a = 2 in Eq. (18) 
according to the classical central limit theorem. This implies -P- S{0 — (O)) and 

consequently 

f(0) = 6(0-{0)). ( 20 ) 


Usually for a hnite system the time a particle spends in a given cell cannot be arbitrarily 
large, and the distribution of sojourn times must be truncated at some value of time. For the 
Hamiltonian mean held model the statistics of sojourn times exhibits a truncated algebraic 
tail with a truncation time diverging with N, implying that the system is non-ergodic only 
in the limit N ^ oo [T8] . 


C. Dynamical functional method 

For a stationary inhnitely divisible processes Y{n), with n integer, the dynamical func¬ 
tional is given by the Fourier transform with unit wave number of the process Y (n) — 
F(0) [35H37|: 

D(n) = (21) 
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where {■ ■ ■) denotes the ensemble average, i. e. an average over many realizations of the 
stochastic process Y{n). The system is ergodic if and only if 

n—1 


k=i 


( 22 ) 


and mixing if and only if: 

lim D(n) = (23) 

n^oo V / I \ / I V / 

We identify Y[n) to the momentum or position variables (or any of their components) at 
time nAt for a fixed time interval At, such that 


1 ^ 

^i[y(n)-y(0)]\ _ _ ^ gi[pfe(nAt)-ps,(0)]^ 

^ k=l 


(24) 


and similarly for the position variable. 

It is useful to dehne the new functions: 


B(n)=C(n)-Ke"'>"')r, 


n—1 


Q(«) = - E E(.k). 


n 


k=l 


Then Eqs. (22) and (23) are rewritten, respectively, as: 


lim Q(n) = 0, 


( 26 ) 


( 26 ) 


and 


lim Ein) = 0. (27) 

Testing for mixing is of special interest as many theories for violent relaxation presuppose a 
good mixing in the initial evolution [6] (see also [56] and references therein). 

We dehne now the time for ergodicity as the time scale such that the difference in time 
and ensemble averages are negligible, and can be estimated using the direct method, which 
simply implements the dehnition of ergodicity. The truncation of sojourn times and the value 


of n in Eqs. (26) and (27) at which E{n) and Q{n) are close to zero, although somehow 
related to the time of ergodicity, are not a priori the same. 


III. SHORT-RANGE INTERACTING SYSTEM - HARD-DISCS IN 2D 


As an illustration of the methods described in the previous section and for comparison 
purposes, let us consider a two-dimensional system composed of identical hard-discs of unit 






mass and diameter d with instantaneous (contact) elastic collisions. Simulations are per¬ 
formed using an event-driven algorithm |1D]. The system is prepared in a homogeneous 
distribution in a square box with sides of length L = 1 and a Gaussian (equilibrium) initial 
distribution of velocities with inverse temperature (3. For economy of space we only per¬ 
form for this model ergodicity tests for the dynamics of the x-component of the momentum 
variable, as the position variables have a similar behavior. To apply the direct method we 
compute the time average of the momentum of each particle and then compute the standard 
deviation of these averages over all N particles, as given by Eq. Results for 

N = 900, N = 2500 and N = 10 000, disc diameter d = 0.002 and (3 = 10.0 are shown in 
Fig. [J. The dependence of the tail of (f) on time and N can be determined by consider¬ 
ing that increments momentum due to collisions are uncorrelated, implying an asymptotic 
time dependence in as a consequence of the classical central limit theorem, and that 
the collisions frequency is proportional to N^ i. e that ~ 1/N^/t as is visible in Fig.[^. 

The thermodynamic limit N —)■ cx), —)■ cxd with N/L'^ = constant, can be inferred from 
Figs. and where some simulations with N = 900 and some values of L are shown. It 
becomes clear that crp^{t) ~ L. From the former results we conclude that: 

N 

lim (Jio (t) = 0, — = constant. (28) 


Therefore, the time for ergodicity vanishes in the thermodynamic limit. 

The average thermal velocity can be dehned at equilibrium by tJ = {v^) = l/\/^,and 

two different time-scales can be dehned from it: the average time between two collisions for 
a given particle Tcoi = l/v, with I the mean free path which is roughly proportional to \/iV, 
and the crossing time, i. e. the time for a particle with average thermal velocity to cross 
the box Ter = 1/F. For the state considered in Fig. [^we have v ~ 0.32. The time required 
for (t) to decrease by two orders of magnitude is roughly two order of magnitudes bigger 
than the crossing time Ter = a/IO. 

Sojourn times are obtained by considering a momentum cell size of Ap^ = 10“^. The 
histograms (distribution functions) of sojourn times are shown in Fig. for N = 2500 
and = 10 000. An exponential tail is evident indicating that the system is ergodic (no 
power law in the tail). Finally for the dynamical functional method. Fig. shows the real 
and imaginary parts of E{n) and Q{n). Were we set n = t for integral values of t. The 


asymptotic values rapidly tend to zero in accordance with Eqs. (26) and (27) (up to some 
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FIG. 1. a) Standard deviation cjp in Eq. (11) of the time average of the x-component of the 


momentum of the computed over all particles for a gas of hard-discs for N = 900, N = 2500 and 
N = 10 000, with disc diameter d = 0.002 and inverse temperature /? = 10.0. The dotted line is the 
curve proportional to ll\/t. b) Same as (a) but with time rescaled as t —)■ t x A^/10 000. The good 
data collapse indicates that o'p^{t) ~ c) Standard deviation for N = 900 but different 

values for the length of the side L of the square box with d = 0.002. d) Same as (c) but with time 
rescaled as t —)• t/L. The good data collapse indicates that crp^{t) ~ L. 


fluctuations), indicating that the system is not only ergodic but also mixing. 


IV. THE HAMILTONIAN MEAN-FIELD MODEL REVISITED 

To better grasp how diverse is the behavior of systems with long-range interactions, we 
briefly discuss the ergodicity of the widely studied HMF model. Its ergodic properties were 
studied in Ref. [18] using the sojourn times method only for homogeneous non-equilibrium 
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FIG. 2. a) Log-Linear plot of normalized histogram of sojourn times with cells in momentum space 
of width Apx = 10“^ and N = 2500 and total simulation time t/ = 1000.0. b) Same as (a) but 
with N = 10 000 and t/ = 237.0. Both histograms show clearly an exponential tail. 




FIG. 3. (Color online) Left panel: Real and imaginary part of Q{t) for a gas of hard-discs with 
N = 2500 and N = 10 000 for the same initial conditions as in Fig. Right panel: Evolution of 


E(t). In both cases we set n = t in Eq. (25). 


QSS. Here we discuss its ergodic properties for homogeneous and a non-homogeneous states 
at the thermodynamics equilibrium. We again restrict ourselves to analyze momentum vari¬ 
ables as position variables show a similar behavior, leaving a more complete analysis for 
the two-dimensional self-gravitating system in next section. Figure shows the standard 
deviation ap in Eq. for two different values of the total energy per particle at ther¬ 
modynamic equilibrium. For e = 0.4 the inverse temperature is /3 ~ 3.02 and the system 
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is in a non-homogeneous state, while for e = 0.8 the state is homogeneous with /3 ~ 1.67 
and a non-vanishing mean-held force. The thermal average velocities are v ~ 0.58 and 
V ~ 0.77, respectively. We hrst observe that the time required for the standard deviation dp 
to decrease by two orders of magnitude is roughly six order of magnitude bigger than the 
crossing time ter = “271/v, i. e. the time required for a particles with average thermal velocity 
to make a full turn on the circle. Figures and Ef show that cTp(t) scales proportional 
to the number of particles N. This is somehow expected for the non-homogeneous case as 
the collisional contribution to the kinetic equation is proportional to 1/iV, but not for the 
homogeneous case as in this case the collisional integral is proportional to 1/A^^. In fact it 
was show in Ref. SB that the statistical momenta of the one-particle distribution function, 
and therefore the distribution function itself, scale in time with a factor 1/N for the non- 
homogeneous state and 1/iV^ for the homogeneous case. In the Vlasov limit the time for 
ergodicity diverges with N, at variance with what occurs for the short range system in the 
previous section. 

Figure]^ shows the distribution of sojourn times for the same initial states. A truncated 
algebraic tail is visible in the plot, in accordance with Ref. [18]. Results for the dynamical 
functional method are shown in Figs. [^and[^ and we can conclude that the system is not 
only ergodic but also mixing as long as N is hnite. To stress this point Fig. shows that 
Q{t) and E{t) scale with N, and similarly to the previous paragraph, the conditions for 
ergodicity and for mixing are not satished in the V —)■ cx) limit. 


V. TWO-DIMENSIONAL GRAVITY 

For a self-gravitating V-body system with Hamiltonian in Eq. (|^, we consider an initial 
condition with particles distributed in space uniformly on a circular strip, with inner radius 
Ri = 40.0 and outer radius R2 = 50.0, and a uniform (waterbag) distribution for the 
momenta in the interval —po < P < Po with pq = 0.5. The system then evolves for a 
total time tf = 10 000 such that it settles in a QSS, which is then used as initial condition 
for the different tests. Figure shows the Kurtosis of the momentum distribution, dehned 
by /C = {p^)/ (/C = 3 for a Gaussian distribution), for the values of N used in our 
simulations, showing that the statistical state of the system is stationary, up to fluctuations, 
as required. 
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FIG. 4. Standard deviation ap in Eq. 0 for the momentum time averages for the HMF model for 
a few values of N at thermodynamic equilibrium, a) Total energy per particle e = 0.4 (magnetized 
state), b) Same as (a) but the time rescaled as t —)• t x 5000/A^. c) Total energy per particle e = 0.8 
(homogeneous state), d) Same as (c) but the time rescaled as t t x 5000/iV. 


The standard deviation of the time average of the x-components of momentnm and 
position, as defined by Eqs. ( [II| ) and (12), respectively, are shown in Fig. 10 The average 
thermal velocity is of order n ~ 1, and a crossing time can be defined by the time a particle 
with the thermal velocity takes to cross the spatial extent Lqss of the core-halo QSS structnre 
formed after the initial violent relaxation, and shown in Fig. i). By inspection, we obtain 
that Lqss ~ 200. This characteristic length does not change significantly with N and time. 
The crossing time case is then given by ter = Lqss/v ~ 200. Comparing Fig. 10 with Figs. 
and we observe that the time for and to decrease by two orders of magnitude are 
roughly the same, and much smaller than for the HMF model in units of crossing time, but 
still much larger than for the short-range system. 
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FIG. 5. Left Panel: Frequency histogram of sojourn times in momentum space for the HMF model 
for the same homogeneous state as in Fig. with cell width Ap = 0.1, total simulation time 
etf = 10®. Right Panel: histogram for the non-homogeneous state in Fig. |^, cell width Ap = 0.4 
and total simulation time tj = 10^. A truncated algebraic tail is visible in both cases. 


The plots of {t) and cr^(t) coincide after a small transient for all valnes of N. The same 
behavior is observed for one-dimensional self-gravitating systems (not shown here). This 
seems to imply that the time for ergodicity is the same for all N and that two-dimensional 
self-gravitating systems are ergodic even in the N ^ oo limit. A few comments are in need 
on this point. The power law decay of the dispersion of can be explained by the fact 
that the final dispersion of the traveled distance in the x direction — Xk{S^)Y/N 

is ronghly independent of N^ and constant in time, which implies that the exponent p 


in Eq. (15) vanishes after the violent relaxation. Eqnation (16) then implies that Up^ is 
proportional to as shown in Fig. 10 This is at variance to what was obtained for the 


HMF model where diffnsion never satnrates 


There is no similar reasoning for the 


behavior of as already explained in Sec. IIA To explain the independence of Up^ (f) on N 
we note that its time evolntion is dne to flnctnations of the force (for finite N) [57H59] . The 


statistics of the force flnctnations for a three-dimensional self-gravitating system were hrst 
obtained by Chandrasekhar [SS] and for one and two-dimensional gravity by Chavanis [6U] . 
Chavanis showed that the flnctnation of the force scale with the density of mass, which is 
proportional to N. Since he does not nse a Kac factor in his expressions, we conclnde for the 
potential in Eq. ([^ the flnctnations of the force are independent of N, and thns the same 
being trne for In onr case, if a Kac factor is not introdnced, then the observation 
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FIG. 6. Dynamical functional method for the HMF model in a homogeneous equilibrium state for 
the same state as in Fig. |^. a) Re[(3(t)]; b) Im[Q(t)]; c) Re[£'(t)]; d) Im[Fi(t)]. In all cases we 
used n = t with integer values of t. 


time required for the system be ergodic diverges with y/N (see discussion after Eq. i. e. 
the system is strictly non-ergodic only in the N ^ oo limit. With a Kac factor the system is 
always ergodic and mixing, although a considerable observation time is required. It is also 
worth noticing that the scaling of cxp^ with N is not the same as for the dynamics of the 
one-particle distribution function, which scales with N. 

Sojourn times are again obtained by diving the momentum space into equal size cells 
of hxed width, associated with discrete states. First we consider cells in momentum and 


position X components, with width Ax = 0.2 and Ap^ = 0.2, respectively. Figure H shows 
the frequency histograms for sojourn times G, which is truncated at quite a small value. 
This truncation is expected to occur for hnite systems as discussed in Ref. [18]. On the 
other hand the truncation values does not vary signihcantly with N, at variance with what 
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FIG. 7. Dynamical functional method for the HMF model in a non-homogeneous equilibrium state 
for the same state as in Fig. |^. a) Re [(5(t)]; b) Im [(5(t)]; c) Re d) Im [E{t)]. 


was observed for the HMF model. An algebraic tail is roughly visible for the intermediate 
values of tg, and deviates from it near the truncation of the distribution, were the statistics 
is less significant. 


Results for the dynamical functional method are given in Figures 12 and showing 


Q{n) and E{n), with n = f for integer values of t. In all cases Eqs. (26) and (27) are 
satished only after a very long time (compare with Fig. which is roughly the same for all 
values of N, in agreement with our results above. 


VI. CONCLUDING REMARKS 

We discussed the ergodic and mixing properties of a two-dimensional self-gravitating 
system using the methods described in Section [TTj The potential of the system is regularized 
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FIG. 8. Left and right panels are the same as in Figs. [6}i and respectively, but with time 
rescaled as t —)> t x iV/5000 in both cases. For the homogeneous case (Left panel) the data collapse 
is very good,while for the inhomogeneous case (right panel), data collapse is only approximate. 



FIG. 9. a) Kurtosis for the momentum distribution as a function of time for the values of N used 
below, b) Positions of particles in the configuration used as initial conditions for simulations for 
N = 4096. 


by the Kac factor in order to have a well dehned thermodynamic (or Vlasov) limit. In order 
to better understand our results we also applied the same tests to a gas of two-dimensional 
hard-discs, and to the HMF model. We showed that all these models are ergodic and mixing, 
although for the HMF model this is true only if N is hnite, as the observation time required 
for the system to be ergodic (and mixing) diverges with N. Nevertheless, for two-dimensional 
gravity, this time scale seems to be independent on N. This is related to the conhning nature 
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FIG. 10. (color online) a) Standard deviation in Eq. (11) for the variables px of each particle for the 


two-dimensional self-gravitating system, integration time step At = 0.05 and softening parameter 
e = 10“^. b) Standard deviation of x. The least-squares fit for the algebraic tail of crp^(t) is given 
by 48.4 X and 2271 x t“0.99 ax{t). 



FIG. 11. a) Sojourn times for cells in momentum x component of momentum space with width 
Apx = 0.2 for some different values of N and total simulation time t/ = 10^. b) Same as (a) but 
with cells in the x component of the position with width Ax = 0.2 and tf = 10^. 

of the interaction potential and the scaling of fluctuations of the force around its mean-held 
value on N. If a Kac factor is not introduced in the potential, then the time for ergodicity 
diverges with increasing N. For the hard-disc gas, the time for ergodicity goes to zero in the 
thermodynamic limit, at variance with both long-range interacting systems studied here. 
The direct and dynamical functional methods proved to be relevant tools for testing 
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FIG. 12. a) Real and b) imaginary part of Q(n) in Eq. (26) for the same simulation and number 


of particles N as in Fig. 11 for the x component of the momentum as stochastic variable, c) Real 


and d) imaginary part of E{n) in Eq. (27). 


ergodicity for the self-gravitating system, while the sojourn times approach gave more am¬ 
biguous results, as the tail of their distribution is too short in the present case to asses if 
there is an algebraic part. 

We note that up to the authors knowledge, the only previous study where such character¬ 
istic time was estimated for long-range interacting systems were devoted to the Hamiltonian 
mean held model, where a N dependence was reported [131 HH] • This result is reobtained 
here for an equilibrium initial condition. 

We also point out that a self-gravitating system with —)■ oo is weakly non-ergodic since 
in the one-particle phase space the constant energy hypersurface is always connected, and 
all regions are necessarily accessible by the dynamics. In this way, a mean-held system can 
strongly non-ergodic if the corresponding energy hypersurface can be composed, at least for 
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FIG. 13. Same as Fig. 12 but For the x component of position as stochastic variable. 


some energy values, by disjoint sets. 
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